%% fig5b.m
% 用于绘制分岔图和最大李雅普诺夫指数图（参考论文 Fig.15a）
% 主程序调用封装的 cnn.m 与 LE.m 函数

clear; clc; close all;

%% 参数设置
a     = 0.7;
alpha = 2.1;      
beta  = 2;
gamma = 2;
Im    = 0.15;        % 外激励幅值
F     = 0.15;        % 外激励频率

% 控制参数 a 的取值范围
b_vals = linspace(0.35, 0.65, 200);

% 时间设置
T_final = 1200;   % 总积分时间
T_trans = 200;   % 舍去瞬态时间

% 分岔图用 ODE 选项
optsODE = odeset('RelTol',1e-9,'AbsTol',1e-9);

% 分段积分参数用于 LE 计算
segment = 5;      % 每段归一化时间
% 总 LE 积分段数由 LE.m 内部计算

% 预分配存储数组
bif_data = cell(length(b_vals),1);  % 用于存储各 a 下分岔图数据（散点）
LE_vals  = zeros(length(b_vals),1); % 存储每个 a 对应的最大李雅普诺夫指数

% 对于每个 a 值计算分岔图数据及 LE 值
parfor i = 1:length(b_vals)
    b = b_vals(i);
    
    % 初始状态：主系统 Y0 = [0;0;0]，扰动向量取非零向量归一化
    Y0 = [0; 0; 0];
    v0 = [1; 0; 0]; v0 = v0/norm(v0);
    X0 = [Y0; v0];
    
    % 1. 分岔图数据计算：积分 [T_trans, T_final]（只提取主系统部分）
    tspan_bif = [T_trans, T_final];
    [t_bif, X_bif] = ode45(@(t,X) cnn(t, X, a, b, alpha, beta, gamma, Im, F), tspan_bif, X0, optsODE);
    Y_bif = X_bif(:,1:3);
    q     = Y_bif(:,1);
    sigma = Y_bif(:,2);
    % 计算 V_x = (alpha + beta*sigma).*q
    Vx = (alpha + beta.*sigma).*q;
    % 利用 findpeaks 提取局部峰值
    [pks, locs] = findpeaks(Vx, t_bif, ...
    'MinPeakProminence', 0.01, ...  % 调整为适当数值
    'MinPeakDistance', 5, ...       % 视采样频率/时间步
    'MinPeakHeight', 0.05);         % 视波形幅值大小
    b_data{i}.b = b * ones(length(pks),1);
    bif_data{i}.Vx  = pks;
    
    % 2. 调用封装函数 LE.m 计算最大李雅普诺夫指数
    LE_vals(i) = LE(a, b, alpha, beta, gamma, Im, F, T_trans, T_final, segment);
end

% 合并所有 a 值下的分岔图数据
all_b_bif   = [];
all_Vx_peaks = [];
for i = 1:length(b_vals)
    all_b_bif   = [all_b_bif; b_data{i}.b];
    all_Vx_peaks = [all_Vx_peaks; bif_data{i}.Vx];
end

%% 绘图
figure;

% 子图1：分岔图
subplot(2,1,1);
scatter(all_b_bif, all_Vx_peaks, 10, 'filled');
xlabel('控制参数 a');
ylabel('V_x 局部峰值');
title('分岔图');
grid on;
xlim([min(b_vals) max(b_vals)]);

% 子图2：最大李雅普诺夫指数图
subplot(2,1,2);
plot(b_vals, LE_vals, 'b.-','LineWidth',1.5);
xlabel('控制参数 a');
ylabel('最大李雅普诺夫指数');
title('李雅普诺夫指数');
grid on;
xlim([min(b_vals) max(b_vals)]);
